// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

#include <vtkSphericalTransform.h>
#include <vtkGeoProjection.h>
#include <vtkGeoTransform.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkStructuredGrid.h>
#include <vtkRectilinearGrid.h>
#include <vtkVisItUtility.h>
#include <DebugStream.h>
#include <vtkCellArrayIterator.h>


// ************************************************************************* //
//  File: avtCartographicProjectionFilter.C
// ************************************************************************* //

#include <avtCartographicProjectionFilter.h>


// ****************************************************************************
//  Method: avtCartographicProjectionFilter constructor
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

avtCartographicProjectionFilter::avtCartographicProjectionFilter()
{
  //this->ProjType = CartographicProjectionAttributes::Hammer;
}


// ****************************************************************************
//  Method: avtCartographicProjectionFilter destructor
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
//  Modifications:
//
// ****************************************************************************

avtCartographicProjectionFilter::~avtCartographicProjectionFilter()
{
}


// ****************************************************************************
//  Method:  avtCartographicProjectionFilter::Create
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

avtFilter *
avtCartographicProjectionFilter::Create()
{
    return new avtCartographicProjectionFilter();
}


// ****************************************************************************
//  Method:      avtCartographicProjectionFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

void
avtCartographicProjectionFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const CartographicProjectionAttributes*)a;
}


// ****************************************************************************
//  Method: avtCartographicProjectionFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtCartographicProjectionFilter with the given
//      parameters would result in an equivalent avtCartographicProjectionFilter.
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

bool
avtCartographicProjectionFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(CartographicProjectionAttributes*)a);
}


// ****************************************************************************
//  Method: avtCartographicProjectionFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the CartographicProjection filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
//  Modifications:
//    Kathleen Biagas, Thu Oct 18 10:09:56 PDT 2012
//    Preserve coordinate data type.
//
//    Jean Favre, Mon Jan  7 17:02:19 CET 2013
//    Added a Central Meridian
//    Added support for PolyData (typically, they are shapefiles for continent lines)
//
//    Eric Brugger, Wed Jul 23 08:30:25 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Kathleen Biagas, Wed Nov  8 10:12:54 PST 2023
//    Fix a couple of issue with VTK 9: numPts should be vtkIdType, not int.
//    Use CellArrayIterator with VTK 9 intead of GetNextCell.
//
// ****************************************************************************

avtDataRepresentation *
avtCartographicProjectionFilter::ExecuteData(avtDataRepresentation *in_dr)
{
/* The projection mappings in VTK, based on the libproj4 library, expect data with coordinates
*  in the range [-180, 180] for longitude, and [-90, 90] for latitude.
*  Further, mapped grids cannot sit on Rectilinear Grids. So, we start by converting rectilinear grids
*  to structured grids (i.e. fully curvilinear)
*  A new grid is created based on the input, coordinates are projected, then remapped to the classic range
*  of [-180, 180] for longitude, and [-90, 90] for latitude
*/
  //
  // Get the VTK data set.
  //
  vtkDataSet *in_ds = in_dr->GetDataVTK();

  int do_type = in_ds->GetDataObjectType();
  int dims[3];
  vtkIdType numPts = in_ds->GetNumberOfPoints();
  vtkPointSet *ds;
  double in_bounds[6], out_bounds[6], tol = 10., p_pt[3], c_pt[3];

  vtkPoints *inPts = vtkVisItUtility::NewPoints(in_ds);
  inPts->Allocate(numPts);
  inPts->SetNumberOfPoints(numPts);

  for(int i=0; i < numPts; i++)
    {
    in_ds->GetPoint(i, c_pt);
// the next two lines used for our Geo-physics grid
    //c_pt[0] = (c_pt[2] - M_PI)*179.999/M_PI;
    //c_pt[1] = (c_pt[1] - M_PI_2)*89.999/M_PI_2;
// this is used for our Geo-physics grid
    c_pt[2] = 0.0;
    inPts->SetPoint(i, c_pt);
    }

  inPts->GetBounds(in_bounds);
  debug4 << "Input Bounds[6] = " << in_bounds[0] << " "  << in_bounds[1] << " " << in_bounds[2] << " " << in_bounds[3] << " " << in_bounds[4] << " " << in_bounds[5] <<  endl;

  vtkGeoProjection *proj = vtkGeoProjection::New();
  proj->SetCentralMeridian(atts.GetCentralMeridian());
  debug4 << "Central Meridian = " << proj->GetCentralMeridian() <<   endl;
  vtkGeoTransform *geoXform = vtkGeoTransform::New();

// Here we could have as many of the 100+ projections available.
// I'm not sure people want to see a pull-down list with 100+ items

  proj->SetName(atts.ProjectionID_ToString(atts.GetProjectionID()).c_str());
  debug4 << "setting name of projection " << atts.ProjectionID_ToString(atts.GetProjectionID()) <<   endl;
  geoXform->SetDestinationProjection(proj);

  vtkPoints *newPoints = vtkPoints::New(inPts->GetDataType());
  geoXform->TransformPoints(inPts, newPoints);

  newPoints->GetBounds(out_bounds);
  debug4 << "Output Bounds[6] = " << out_bounds[0] << " "  << out_bounds[1] << " " << out_bounds[2] << " " << out_bounds[3] << " " << out_bounds[4] << " " << out_bounds[5] <<  endl;

// rescale coordinates to the original set of longitude and latitude ranges.
// this is motivated by the fact that the continents shapefile provided by Brad for the SC'12 tutorial
// does not range fully to the North Pole. It was wrong to stretch its projection to +90 degrees.
  double alpha;
  for(int i=0; i < numPts; i++)
    {
    newPoints->GetPoint(i, c_pt);
    alpha = (c_pt[0] - out_bounds[0]) / (out_bounds[1]-out_bounds[0]);
    c_pt[0] = alpha * in_bounds[1] + (1.0 - alpha) * in_bounds[0];
    alpha = (c_pt[1] - out_bounds[2]) / (out_bounds[3]-out_bounds[2]);
    c_pt[1] = alpha * in_bounds[3] + (1.0 - alpha) * in_bounds[2];
    // only touch the X and Y coordinates. Ignore Z
    newPoints->SetPoint(i, c_pt);
    }
  // mark the coordinates as Modified to force a call to vtkPoints::ComputeBounds()
  // otherwise a ResetView() has trouble resetting.
  newPoints->Modified();
  newPoints->GetBounds(out_bounds);
  debug4 << "Output Bounds[6] = " << out_bounds[0] << " "  << out_bounds[1] << " " << out_bounds[2] << " " << out_bounds[3] << " " << out_bounds[4] << " " << out_bounds[5] <<  endl;

  vtkCellArray *ca_n;
  int changeOfSigns = 0, k;

  switch(do_type) {
    case VTK_STRUCTURED_GRID:
      ds = vtkStructuredGrid::New();
      static_cast<vtkStructuredGrid *>(in_ds)->GetDimensions(dims);
// make sure 3rd dimension is set in case we were given 2D datasets (e.g pressure.cdf)
      dims[2] = 1;
      static_cast<vtkStructuredGrid *>(ds)->SetDimensions(dims);
      ds->ShallowCopy(in_ds);
    break;
    case VTK_RECTILINEAR_GRID:
      ds = vtkStructuredGrid::New();
      static_cast<vtkRectilinearGrid *>(in_ds)->GetDimensions(dims);
// make sure 3rd dimension is set in case we were given 2D datasets (e.g pressure.cdf)
      dims[2] = 1;
      static_cast<vtkStructuredGrid *>(ds)->SetDimensions(dims);
      ds->ShallowCopy(in_ds);
    break;
    case VTK_POLY_DATA:
    { // new scope
      // some special treatment is done here for polylines which -
      // when projected -  "fall on the other side of the Earth".
      // Detect an line segment within the polyline which has a very
      // long length and split.
      ds = vtkPolyData::New();
      ca_n = vtkCellArray::New();
      vtkPolyData *pd = vtkPolyData::SafeDownCast(ds);
      pd->SetPolys(ca_n);
      ca_n->Delete();

      vtkIdType npts;

      const vtkIdType *pts=nullptr;

      auto ca = vtk::TakeSmartPointer(pd->GetPolys()->NewIterator());
      for (ca->GoToFirstCell(); !ca->IsDoneWithTraversal(); ca->GoToNextCell())
      {
          ca->GetCurrentCell(numPts, pts);

         // for each polygon, change for big changes in coordinates and split lines
         changeOfSigns = 0;
         k = npts-1;
         // start from end and split if necessary
         for(int j =npts-1; j >0; j--)
           {
           newPoints->GetPoint(pts[j-1], p_pt); // previous pt
           newPoints->GetPoint(pts[j], c_pt);    // current pt
           // compute a 2d distance
           if(sqrt((p_pt[0] - c_pt[0])*(p_pt[0] - c_pt[0]) + (p_pt[1] - c_pt[1])*(p_pt[1] - c_pt[1])) > tol )
             {
             changeOfSigns++;
             ca_n->InsertNextCell(k-j+1, &pts[j]); k = j-1;
             }
           }
        if(changeOfSigns == 0)
          {
          ca_n->InsertNextCell(npts, pts);
          } // full polygon (original)
        else
          {
          ca_n->InsertNextCell(k+1, pts); // what is left-over after all the splits
          }
         }
    }
    break;
    default:
      debug4 << "not supported for this grid type"  <<endl;
      return in_dr;
  }

  ds->SetPoints(newPoints);

  if(do_type != VTK_POLY_DATA)
  {
    debug4 << "dims = " << dims[0] << " x "  << dims[1] << " x " << dims[2] <<  endl;
  }

  if(do_type == VTK_RECTILINEAR_GRID)
    inPts->Delete();

  newPoints->Delete();
  proj->Delete();
  geoXform->Delete();

  avtDataRepresentation *out_dr = new avtDataRepresentation(ds,
      in_dr->GetDomain(), in_dr->GetLabel());

  ds->Delete();

  return out_dr;
}

// ****************************************************************************
// Method: avtCartographicProjectionFilter::UpdateDataObjectInfo
//
// Purpose:
//   Update the data object information.
//
// Programmer: Jean M. Favre
// Creation:   Mon Apr  7 16:11:17 PDT 2014
//
// Modifications:
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtCartographicProjectionFilter::UpdateDataObjectInfo(void)
{
  avtDataAttributes &dataAtts = GetOutput()->GetInfo().GetAttributes();
  dataAtts.SetXUnits("");
  dataAtts.SetXLabel("Longitude");
  dataAtts.SetYUnits("");
  dataAtts.SetYLabel("Latitude");
  dataAtts.SetSpatialDimension(2);

  dataAtts.AddFilterMetaData("CartographicProjection");
}
